3. FEM-BEM Coupling#
The ngbem boundary element addon project initiated by Lucy Weggeler (see https://weggler.github.io/ngbem/intro.html) is now partly integrated into core NGSolve. Find a short and sweet introduction to the boundary element method there.
In this demo we simulate a plate capacitor on an unbounded domain.
from ngsolve import *
from netgen.occ import *
from ngsolve.solvers import GMRes, MinRes
from ngsolve.webgui import Draw
from ngsolve.bem import *
largebox = Box ((-2,-2,-2), (2,2,2) )
eltop = Box ( (-1,-1,0.5), (1,1,1) )
elbot = Box ( (-1,-1,-1), (1,1,-0.5))
largebox.faces.name = "outer" # coupling boundary
eltop.faces.name = "el1" # Dirichlet boundary
elbot.faces.name = "el2" # Dirichlet boundary
eltop.edges.hpref = 1
elbot.edges.hpref = 1
shell = largebox-eltop-elbot # FEM domain
shell.solids.name = "air"
mesh = shell.GenerateMesh(maxh=0.8)
mesh.RefineHP(2)
ea = { "euler_angles" : (-67, 0, 110) }
Draw (mesh, clipping={"x":1, "y":0, "z":0, "dist" : 1.1}, **ea);
On the exterior domain \(\Omega^c\), the solution can be expressed by the representation formula:
\[ x \in \Omega^c: \quad u(x) = - \int\limits_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{1}{\| x-y\|} } \, \gamma_1 u (y)\, \mathrm{d}\sigma_y + \int\limits_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{\langle n(y) , x-y\rangle }{\| x-y\|^3} } \, \gamma_0 u (y)\, \mathrm{d}\sigma_y\,, \]
where \(\gamma_0 u = u\) and \(\gamma_1 u = \frac{\partial u}{\partial n}\) are Dirichlet and Neumann traces. These traces are related by the Calderon projector
\[ \begin{align}\begin{aligned}\begin{split} \left( \begin{array}{c} \gamma_0 u \\ \gamma_1 u \end{array}\right) = \left( \begin{array}{cc} -V & \frac12 + K \\ \frac12 - K^\intercal & -D \end{array} \right) \left( \begin{array}{c} \gamma_1 u \\ \gamma_0 u \end{array}\right) $$.\end{split}\\The $V$, $K$ are the single layer and double layer potential operators, and $D$ is the hypersingular operator.\\On the FEM domain we have the variational formulation\end{aligned}\end{align} \]
\int_{\Omega_\text{FEM}} \nabla u \nabla v , dx - \int_\Gamma \gamma_1 u v , ds = 0 \qquad \forall , v \in H^1(\Omega_\text{FEM}) $$
We use Calderon’s represenataion formula for the Neumann trace:
\[
\int_{\Omega_\text{FEM}} \nabla u \nabla v \, dx - \int_\Gamma \left( \left( \tfrac{1}{2} - K^\intercal\right) \,\gamma_1 u - D \, \gamma_0 u\right) v = 0 \qquad \forall \, v \in H^1(\Omega_\text{FEM})
\]
To get a closed system, we use also the first equation of the Calderon equations. To see the structure of the discretized system, the dofs are split into degrees of freedom inside \(\Omega\), and those on the boundary \(\Gamma\). The FEM matrix \(A\) is split accordingly. We see, the coupled system is symmetric, but indefinite:
\[\begin{split} \left( \begin{array}{ccc } A_{\Omega\Omega} & A_{\Omega\Gamma} & 0 \\ A_{\Gamma\Omega} & A_{\Gamma\Gamma } + D & -\frac12 M^\intercal + K^\intercal \\ 0 & -\frac12 M + K & -V \end{array}\right) \left( \begin{array}{c} u \\ \gamma_0 u \\ \gamma_1 u \end{array}\right) = \left( \begin{array}{c} F_{\Omega} \\ F_{\Gamma}\\ 0 \end{array}\right) \,. \end{split}\]
Generate the finite element space for \(H^1(\Omega)\) and set the given Dirichlet boundary conditions on the surfaces of the plates:
order = 4
fesH1 = H1(mesh, order=order, dirichlet="el1|el2")
print ("H1-ndof = ", fesH1.ndof)
H1-ndof = 91269
The finite element space \(\verb-fesH1-\) provides \(H^{\frac12}(\Gamma)\) conforming element to discretize the Dirichlet trace on the coupling boundary \(\Gamma\). However we still need \(H^{-\frac12}(\Gamma)\) conforming elements to discretize the Neumann trace of \(u\) on the coupling boundary. Here it is:
fesL2 = SurfaceL2(mesh, order=order-1, dual_mapping=True, definedon=mesh.Boundaries("outer"))
print ("L2-ndof = ", fesL2.ndof)
L2-ndof = 4035
fes = fesH1 * fesL2
u,dudn = fes.TrialFunction()
v,dvdn = fes.TestFunction()
a = BilinearForm(grad(u)*grad(v)*dx, check_unused=False).Assemble()
gfudir = GridFunction(fes)
gfudir.components[0].Set ( mesh.BoundaryCF( { "el1" : 1, "el2" : -1 }), BND)
f = LinearForm(fes).Assemble()
res = (f.vec - a.mat * gfudir.vec).Evaluate()
Generate the the single layer potential \(V\), double layer potential \(K\) and hypersingular operator \(D\):
n = specialcf.normal(3)
with TaskManager():
V = LaplaceSL(dudn*ds("outer"))*dvdn*ds("outer")
K = LaplaceDL(u*ds("outer"))*dvdn*ds("outer")
D = LaplaceSL(Cross(grad(u).Trace(),n)*ds("outer"))*Cross(grad(v).Trace(),n)*ds("outer")
M = BilinearForm(u*dvdn*ds("outer"), check_unused=False).Assemble()
Setup the coupled system matrix and the right hand side:
sym = a.mat+D.mat - (0.5*M.mat+K.mat).T - (0.5*M.mat+K.mat) - V.mat
rhs = res
bfpre = BilinearForm(grad(u)*grad(v)*dx+1e-10*u*v*dx + dudn*dvdn*ds("outer") ).Assemble()
pre = bfpre.mat.Inverse(freedofs=fes.FreeDofs(), inverse="sparsecholesky")
Compute the solution of the coupled system:
with TaskManager():
sol_sym = MinRes(mat=sym, rhs=rhs, pre=pre, tol=1e-8, maxsteps=200, printrates=True)
LinearSolver iteration 1, residual = 18.46865730865644
LinearSolver iteration 2, residual = 2.340062830525568
LinearSolver iteration 3, residual = 2.082508582473063
LinearSolver iteration 4, residual = 0.42147022859368977
LinearSolver iteration 5, residual = 0.26625863647556525
LinearSolver iteration 6, residual = 0.11205519202260256
LinearSolver iteration 7, residual = 0.057861992879118244
LinearSolver iteration 8, residual = 0.035293799810918056
LinearSolver iteration 9, residual = 0.015622853897262664
LinearSolver iteration 10, residual = 0.014077891549896164
LinearSolver iteration 11, residual = 0.014053434842592695
LinearSolver iteration 12, residual = 0.003550247793548699
LinearSolver iteration 13, residual = 0.003096966218249709
LinearSolver iteration 14, residual = 0.002083724623773071
LinearSolver iteration 15, residual = 0.0016137942053991209
LinearSolver iteration 16, residual = 0.0013415743220941197
LinearSolver iteration 17, residual = 0.0013275776848339148
LinearSolver iteration 18, residual = 0.0007480874464816323
LinearSolver iteration 19, residual = 0.0007138749384091383
LinearSolver iteration 20, residual = 0.00057019812483812
LinearSolver iteration 21, residual = 0.000563713804681432
LinearSolver iteration 22, residual = 0.000383788901185162
LinearSolver iteration 23, residual = 0.00038346135532268206
LinearSolver iteration 24, residual = 0.0002819115766637901
LinearSolver iteration 25, residual = 0.0002735291921814488
LinearSolver iteration 26, residual = 0.0002151366536708616
LinearSolver iteration 27, residual = 0.00020073895249888356
LinearSolver iteration 28, residual = 0.00013795110697862817
LinearSolver iteration 29, residual = 0.00013767309251228196
LinearSolver iteration 30, residual = 0.00012064207804577526
LinearSolver iteration 31, residual = 0.00010761004002411436
LinearSolver iteration 32, residual = 7.815669958088078e-05
LinearSolver iteration 33, residual = 7.766309129369102e-05
LinearSolver iteration 34, residual = 6.963167856977043e-05
LinearSolver iteration 35, residual = 6.197479416126628e-05
LinearSolver iteration 36, residual = 6.110378466223794e-05
LinearSolver iteration 37, residual = 4.408912650385276e-05
LinearSolver iteration 38, residual = 4.003428990753054e-05
LinearSolver iteration 39, residual = 3.706913827958613e-05
LinearSolver iteration 40, residual = 3.060671629523218e-05
LinearSolver iteration 41, residual = 2.7262324111116e-05
LinearSolver iteration 42, residual = 2.311872091682488e-05
LinearSolver iteration 43, residual = 2.2187834117085572e-05
LinearSolver iteration 44, residual = 2.1062267392025543e-05
LinearSolver iteration 45, residual = 1.9111099817039495e-05
LinearSolver iteration 46, residual = 1.9075210465238373e-05
LinearSolver iteration 47, residual = 1.74977232406471e-05
LinearSolver iteration 48, residual = 1.574879534670412e-05
LinearSolver iteration 49, residual = 1.2286066130723375e-05
LinearSolver iteration 50, residual = 1.2239528725466086e-05
LinearSolver iteration 51, residual = 1.0054445260235162e-05
LinearSolver iteration 52, residual = 1.0036035629501306e-05
LinearSolver iteration 53, residual = 7.963163124668724e-06
LinearSolver iteration 54, residual = 7.96271087227098e-06
LinearSolver iteration 55, residual = 7.694618539792837e-06
LinearSolver iteration 56, residual = 6.7535722988566266e-06
LinearSolver iteration 57, residual = 6.560603126335692e-06
LinearSolver iteration 58, residual = 5.776177872134021e-06
LinearSolver iteration 59, residual = 5.683864380734816e-06
LinearSolver iteration 60, residual = 4.676861662449405e-06
LinearSolver iteration 61, residual = 4.606728802970939e-06
LinearSolver iteration 62, residual = 3.963055733921582e-06
LinearSolver iteration 63, residual = 3.660259877144058e-06
LinearSolver iteration 64, residual = 3.4623689512087837e-06
LinearSolver iteration 65, residual = 3.393967594604926e-06
LinearSolver iteration 66, residual = 2.7626999825999792e-06
LinearSolver iteration 67, residual = 2.76267358440659e-06
LinearSolver iteration 68, residual = 2.461990411606185e-06
LinearSolver iteration 69, residual = 2.4494580218698885e-06
LinearSolver iteration 70, residual = 2.1620228169065188e-06
LinearSolver iteration 71, residual = 2.0947239732977038e-06
LinearSolver iteration 72, residual = 1.7296792574154861e-06
LinearSolver iteration 73, residual = 1.6832748746496014e-06
LinearSolver iteration 74, residual = 1.6219706651074022e-06
LinearSolver iteration 75, residual = 1.4388437428061557e-06
LinearSolver iteration 76, residual = 1.4388431338620584e-06
LinearSolver iteration 77, residual = 1.1051343600676078e-06
LinearSolver iteration 78, residual = 1.0320274044996104e-06
LinearSolver iteration 79, residual = 9.47117957094968e-07
LinearSolver iteration 80, residual = 9.200861941113419e-07
LinearSolver iteration 81, residual = 8.104438883485559e-07
LinearSolver iteration 82, residual = 8.05228777002287e-07
LinearSolver iteration 83, residual = 7.39581829396173e-07
LinearSolver iteration 84, residual = 6.778707159825941e-07
LinearSolver iteration 85, residual = 5.897108430277277e-07
LinearSolver iteration 86, residual = 5.534413854493185e-07
LinearSolver iteration 87, residual = 5.040212322886497e-07
LinearSolver iteration 88, residual = 4.6821161084395414e-07
LinearSolver iteration 89, residual = 4.4368697035005224e-07
LinearSolver iteration 90, residual = 3.975283306395698e-07
LinearSolver iteration 91, residual = 3.3775057338079296e-07
LinearSolver iteration 92, residual = 3.2299945870887103e-07
LinearSolver iteration 93, residual = 3.2256873811828665e-07
LinearSolver iteration 94, residual = 2.7928663204101156e-07
LinearSolver iteration 95, residual = 2.7644719603003633e-07
LinearSolver iteration 96, residual = 2.621210431213476e-07
LinearSolver iteration 97, residual = 2.4553371791659154e-07
LinearSolver iteration 98, residual = 2.1394705180034105e-07
LinearSolver iteration 99, residual = 1.9893719713486604e-07
LinearSolver iteration 100, residual = 1.716608908084191e-07
gfu = GridFunction(fes)
gfu.vec[:] = sol_sym + gfudir.vec
Draw(gfu.components[0], clipping={"x" : 1, "y":0, "z":0, "dist":0.0, "function" : True }, **ea, order=2);
The Neumann data:
Draw (gfu.components[1], **ea, order=3);